Noise-induced perturbations of dispersion-managed solitons 
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We study noise-induced perturbations of dispersion-managed solitons by developing soliton perturbation theory 
for the dispersion-managed nonlinear Schrodinger (DMNLS) equation, which governs the long-term behavior 
of optical fiber transmission systems and certain kinds of femtosecond lasers. We show that the eigenmodes and 
generalized eigenmodes of the linearized DMNLS equation around traveling-wave solutions can be generated 
from the invariances of the DMNLS equations, we quantify the perturbation-induced parameter changes of the 
solution in terms of the eigenmodes and the adjoint eigenmodes, and we obtain evolution equations for the 
solution parameters. We then apply these results to guide importance-sampled Monte-Carlo simulations and 
reconstruct the probability density functions of the solution parameters under the effect of noise. 

PACS numbers: 42.65.Sf, 42.65.Tg 



I. INTRODUCTION 

Dispersion management has become an essential compo- 
nent not only of modern optical fiber communication sys- 
tems [6, 22], but also of certain femtosecond lasers [37]. The 
performance of both kinds of systems is affected by noise, 
which is an essential source of system failures. Because these 
systems are designed to operate with extremely high accura- 
cies (typical values are one error per 10 12 bits in communi- 
cations and 1 part in 10 18 for lasers used in optical atomic 
clocks), calculating failure rates analytically is extremely dif- 
ficult since failures result from the occurrence of unusually 
large (and therefore atypical) deviations. At the same time, 
direct Monte-Carlo computations of failure rates are impracti- 
cal due to the exceedingly large number of samples that would 
be necessary to obtain a reliable estimate. 

The effect of noise on optical transmission systems modeled 
by the nonlinear Schrodinger (NLS) equation has recently 
been studied [23, 25, 29, 35] using a variance reduction tech- 
nique called importance sampling (IS). In brief, IS biases the 
Monte-Carlo simulations in such a way as to artificially in- 
crease the probability of achieving the rare events of interest, 
while correcting for the bias using appropriate likelihood ra- 
tios (e.g., see Refs. [10, 30]). Use of IS makes it possible to 
efficiently estimate extremely small probabilities. 

The key in successfully applying importance sampling lies in 
biasing towards the most likely noise realizations that lead to 
system failures. In the above-cited works, this was achieved 
by taking advantage of well-known results about the behavior 
of solutions of the NLS equation linearized around a soliton 
solution. This knowledge is not available, however, in sys- 
tems with dispersion management. The aim of this work is to 
address and overcome this problem. We do so by employing 
the dispersion-managed NLS (DMNLS) equation which gov- 
erns the long-term dynamics of dispersion-managed optical 
systems [1, 5, 12]. 

The layout of this paper is as follows. In section II we develop 
a perturbation theory for dispersion-managed systems. First 
we study the connection between the invariances of DMNLS 
equation and solutions of the linearized DMNLS equation. 
We then show how the equation invariances are connected 



to the existence of traveling-wave solutions. We also show 
that, as for the NLS equation, the linearized DMNLS around 
such traveling-wave solution can be expressed in terms of an 
ordinary differential operator whose eigenmodes and general- 
ized eigenmodes can also be generated from the same invari- 
ances. Finally, we use these linear modes and their adjoints to 
quantify the perturbation-induced parameter changes, using 
the relation between the linear modes and the derivatives of 
the solution with respect to the invariance parameters. In sec- 
tion III we use these theoretical results to guide importance- 
sampled Monte-Carlo simulations of noise-induced perturba- 
tions in dispersion-managed lightwave systems and recon- 
struct the probability density functions of the output solution 
parameters. 



II. SYMMETRIES AND PERTURBATIONS OF 
DISPERSION-MANAGED SYSTEMS 



It is well-known that the propagation of coherent optical 
pulses in dispersion-managed systems can be described by the 
following perturbed NLS equation: 



du 



+ \d{t/t a )^-j+g{t/t a )\u\ 2 u = Q. 



(2.1) 



Here all quantities are expressed in dimensionless units; t is 
the propagation distance, x is the retarded time (that is, the 
time in a reference frame that moves with the group veloc- 
ity of the pulse), and u(x,t) is the slowly varying envelope of 
the optical field, rescaled (if necessary in communications) to 
take into account periodic loss and amplification. The func- 
tion d(t I to) represents the local dispersion, while g(t/t a ) de- 
scribes the periodic power variation due to loss and ampli- 
fication. (That is, the optical amplitude is proportional to 
y/g(t/t a ) u(x,t).) Both d( ■ ) and g( ■ ) are taken to be periodic 
with unit period. The particular choice of d(t/t a ) is called 
a dispersion map, and the quantity t a is called the dispersion 
map period. Systems described by Eq. (2.1) include modern 
optical fiber communication systems [6, 13] as well as certain 
femtosecond lasers [37]. 
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A. The DMNLS equation, invariances and soliton solutions 



management, and are given respectively by 



Some of the invariances of the "pure" NLS equation (namely, 
Eq. (2.1) with d(-) = g(-) = 1) are lost with dispersion man- 
agement. (More precisely, time translations, scaling and 
Galilean invariance, although a generalized Galilean invari- 
ance exists.) Moreover, Eq. (2.1) is a nonlinear partial differ- 
ential equation (PDE) with non-constant coefficients which 
contain large and rapid variations; the asymptotic behav- 
ior of its solutions is therefore not apparent. As shown in 
Ref. [1], an appropriate multiple-scale analysis of Eq. (2.1) 
shows that, once the periodic compression-expansion breath- 
ing of the pulse is properly factored out, the core pulse shape 
obeys a nonlinear, nonlocal equation of nonlinear-Schrodinger 
type called the dispersion-managed NLS (DMNLS) equation. 
Without repeating the derivation here, we note that the key 
is to split the dispersion d(t/t a ) into the sum of two compo- 
nents: a mean value d and a term describing the large, zero- 
mean rapid variations corresponding to the large local values 
of dispersion: 



d{t/t a )=d+]-M){t/t a ). 
To leading order, the solution of Eq. (2.1) is then 

u((0,t)=u'{(0,t)e- ic ^ al l 2 , 



where C,~t/t a and 



C(C)=Co + .[AD(CVC, 
o 



(2.2) 



(2.3) 



(2.4) 



with Co an arbitrary integration constant. Above and here- 
after, f((o) = 7 a [/(*)] = / e~ imx f(x) dx is the Fourier trans- 
form of f(x). The exponential factor in Eq. (2.3) accounts 
for the rapid breathing of the pulse, while the slowly varying 
envelope u((0,t) satisfies the DMNLS equation, which in the 
physical and Fourier domains, is respectively (omitting primes 
for simplicity) 



du 

l Tt 



dx 2 



+ II u (x+x') u (x+x") u *x+x'+x") R (x'jc") dx'dx" - 



-to) 



1 1 



(2n) 2 



o 



(2.6a) 



R(x',x") = -JJe- 



r{(o'(o")d(d'd(o" . (2.6b) 



Note that both focusing and defocusing cases can be obtained 
in Eqs. (2.5) depending on the sign of the average disper- 
sion d. 

It is crucial to realize that the DMNLS equation and its solu- 
tions depend implicitly on a parameter, called the dimension- 
less reduced map strength, which quantifies the size of the 
zero-mean dispersion fluctuations. The map strength s can be 
defined for any dispersion map as 



AD|| 1 = if|AD(C)MC- 



(2.7) 



(2.5a) 



One can then formally obtain the dependence of the ker- 
nels r(y) and R(x',x") on the map strength by writing AD(Q 
and C(Q in terms of normalized functions which only de- 
pend on the shape of the zero-mean dispersion variations. 
Namely, one writes AD(Q = 4* ADi (Q and C(%) = 4sC\ (t) , 
where Ci(£) is given by Eq. (2.4) with AD(Q replaced by 
ADi(£). In this way, one can conveniently study cases with 
different map strengths entirely within the framework of the 
DMNLS equation, without needing to refer to Eq. (2.1). Of 
course, in the limit s — ► one obtains r(y) — > l/(27r) 2 , and 
R(x',x") -> d(x')d(x"). That is, as s the DMNLS equa- 
tion (2.5) reduces to the pure NLS equation. 

The DMNLS equation (2.5) is a reduced model that retains the 
essential features of dispersion-managed systems while by- 
passing the complicated dynamics that take place within each 
dispersion map. As such, it has proved to be a useful model 
to investigate the long-time behavior of dispersion-managed 
systems [1-5, 8, 1 1, 12, 17-20, 26, 28, 33, 34]. 

In the case where loss and gain are perfectly balanced (e.g., 
with distributed Raman amplification in communication sys- 
tems [7]) it is g( ■ ) = 1, and both kernels can then be made real 
by proper choice of the integration constants. Here we assume 
that this has been done. Then, in the special but physically im- 
portant case of a piecewise constant, two-step dispersion map, 
the kernels assume a particularly simple form [1]: 



.du ij-2- 

i— idea u 

dt 2 



+ //> 



'(co+(o') u ((o+co") u (co+m'+co") r ( m ' m ") 



dco'dca" = 0, 



(2.5b) 



where the asterisk denotes complex conjugation and where for 
brevity we introduced the shorthand notations U( x ) = u(x,t), 
t5((u) = u(co,t) etc. Throughout this work, integrals are com- 
plete unless limits are explicitly stated. The integration ker- 
nels r(y) and R(x',x") in Eqs (2.5) quantify the average non- 
linearity over a dispersion map mitigated by the dispersion 



1 sinsy 
(2tz;) 2 sy ' 



R(x',x") 



1 



2k\s 



ci(x'x"/s), (2.8) 



where ci(x) = f™ cos(xy) /y dy. Note that in this case both 
kernels are independent of the particular shape of the zero- 
mean dispersion variations. The same kernels, apart from a 
factor 2, also arise for the DMNLS equation as a model of 
certain femtosecond lasers [5]. 

Stationary solutions of the DMNLS equation are obtained by 
looking for solutions of the form: 



u st (x,t;s) =e iX2, l 2 f{x;s). 



(2.9) 
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The Fourier transform /(©) of f(x) then solves the nonlinear 
integral equation [1] 

{X 2 + d(o 2 )f [m) = 

2 JJ f(o}+co')f{a>+m")f(oi+m'+oi") r (co'm") dco'dco" , (2.10) 

which can be efficiently solved numerically (cf. Appendix 
A 2). Note that, like in the NLS equation, the "soliton eigen- 
value" X is also the peak amplitude of the pulse, indepen- 
dently of map strength. 

Remarkably, some of the invariances of the standard NLS 
equation that are destroyed by dispersion management are re- 
covered by the DMNLS equation. In particular, the following 
NLS invariances also hold for the DMNLS equation without 
modification: phase invariance, position invariance, time in- 
variance and Galilean invariance, which are respectively 



u — > e u , 

X ^ X Xq . 

t—>t-t„ 



1-0 J 



x-£lt, 



u 



JSlx-i&t/l 



(2.11a) 
(2.11b) 
(2.11c) 

(2. lid) 



On the other hand, unlike the NLS equation, the DMNLS 
equation is not scale invariant. Rather, the DMNLS equation 
admits a generalized scaling invariance which also involves 
the map strength: 



k/A, 



■t/A 2 



■ S /A 2 , 



u — > Au. 

(2.1 le) 



Starting from Eq. (2.9), a four-parameter family of traveling- 
wave solutions can then be generated by using the invariances 
(2.11): 

u(x,t) = Ae''^ + 3( A2 - a2 >+«/(A(x -x -Clt);A 2 s) , 

(2.12a) 

or, equivalently, 

u(x,t) = Ae i ^ + ^f(A(x-X);A 2 s), (2.12b) 

where 

X(t)=x a + Qt, 4>(f) = (A 2 - £l 2 )t/2 + <p (2.13) 

are respectively the mean position and an overall phase. Note 
that for the stationary solution (2.9), time translations can be 
expressed as a composition of phase transformations and po- 
sition translations. Hence, even though five invariances ex- 
ist, there are only four independent solution parameters for 
traveling-wave solutions. When the kernel r(y) is real, f(x) 
can be taken to be real and even. Solutions (2.12) are usually 
referred as dispersion-managed solitons (DMS). 

B. Linear modes of the DMNLS equation 

We now consider the stability of solutions under perturba- 
tions. If u(x,t) is any solution of the DMNLS equation and 



u(x,t) +ev(x,t) is also a solution, v(x,t) solves the linearized 
DMNLS equation around u(x,t); namely, L[v,u] — 0, where 
[18,26] 

dv ■ -d 2 v 

L\v,u] = -=: jd^r 

11 dt 2 dx 2 

- 2i JJ u (.x+x") u (x+x>+x") v (x+x>) R (x's") dx ' dx " 
~ i fJ u (x+x') U (x+x") V (x+x'+x») R (x's") dx ' dx " ( 2 - 14 > 

is the linearized DMNLS operator. Since the DMNLS equa- 
tion is not integrable [3, 17], its linear modes cannot be 
derived from the inverse scattering method as for the NLS 
equation [15]. The linear modes, however, can be gener- 
ated using the invariances (2.11). Suppose that u(x,t) solves 
some PDE, and consider a generic infinitesimal transforma- 
tion u(x,t) — > u e (x,t), with 



and 



u e (x,t) — u(x,t) + ev(x,t) + 0(e 2 ) , 
du e (x,t) 



v(x,t) = ■ 



de 



(2.15a) 
(2.15b) 



e=0 



If the PDE is invariant under the above transformation, one 
verifies that v(x,t) is in the nullspace of L[- ,«]; namely, v(x,t) 
is a solution of the linearized PDE about the given solu- 
tion. When applied to the DMNLS equation, this construc- 
tion yields four solutions of the linearized DMNLS equation 
around any solution u(x,t). More precisely, these solutions 
of the linearized DMNLS equation, which are associated with 
the phase, distance translation, Galilean invariance and scale 
invariance, are respectively 



V 3 



IXU 



du 
> 

dx ' 



V2 = - 



du 

dx ' 



(2.16a) 



v 4 



du „ du „ du 
u+x— + 2t—+2s—. 
dx dt ds 



(2.16b) 



Note that V3 and V4 are not bounded in time. Using the 
fact that L[v] = w implies L[tv] = tw + v, however, one can 
convert V3 and V4 into bounded elements of the generalized 
nullspace of L. Note also that a further solution of the lin- 
earized DMNLS equation, namely V5 = — u t , can be generated 
from invariance with respect to time translations. This fifth so- 
lution of the linearized equation, however, is not linearly in- 
dependent from the other four if u(x,t) is the traveling-wave 
solution (2.12), since then 



v 5 = \{A 2 + Cl 2 )v l +Q.v 2 . 



(2.17) 



For traveling-wave solutions, it is possible to express the lin- 
earized DMNLS equation in terms of an ordinary differential 
operator by performing a change of coordinates to the comov- 
ing frame (%,t'), with E, = x—X(t) and t' = t , and writing 
u(x,t) and v(x,t) respectively as 



u(x,t)=e l,d U(t;), v(x,t)=e'*y(t;,t>), 



(2.18) 
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FIG. 1 : The shape of the neutral eigenmodes and generalized eigen- 
modes of the DMNLS equation for s = 2 (solid lines), compared to 
the corresponding modes for the NLS equation (s = 0, dashed lines), 
with d = 1 and A = 2 in both cases. Blue lines (top right and bottom 
right plots) show components in phase with the pulse; red lines (top 
left and bottom left plots) show components 90° out of phase. 

where &(x,t) — Qjc + <P(t) is the local phase and U(£) = 
Af(A%) the pulse envelope. Substituting Eq. (2.18) into 
Eq. (2.14) yields 



-i0 



L\v;u] = ^-A\y;U], 



(2.19) 



where 



-d 2 y 



My,u] = - 2 d^f 2 -- 2 Ab 

+ 2i SI u ^H") u kH'H") y ^H') R ^'^") d % d %' 

Following standard terminology, we will call y a neutral eigen- 
mode if A[y] =0 and a generalized eigenmode if A[y] is a 
neutral eigenmode. Using Eqs. (2.18) and (2.19), one can as- 
sociate each solution of the linearized DMNLS in Eqs. (2.16) 
to a neutral eigenmode or generalized eigenmode of A. After 
rearranging terms, one can obtain the following set of modes 
and generalized modes: 



y$, = iU 
1 

ya. = is U , y A = -. 



yx 



du_ 



u 



.du_ 



■2s 



dU 



(2.21a) 
(2.21b) 



(where the subscript associates each mode to the solution pa- 
rameters changed by the transformation), which satisfy the 
following relations: 



Afro] = yx , 



AM 
Abu] = 



■Ay®. 



(2.22a) 
(2.22b) 



Note the explicit dependence on s of the amplitude mode y A , 
as well as, of course, of the corresponding solution of the lin- 
earized DMNLS equation, V4. Note also that, as for the NLS 



equation, some freedom exists in the definition of the gener- 
alized modes, as well as in the normalization of all modes. 
Importantly, for Q = A , X , 4> it is 

du 



dQ 



= e m y Q , 



while 



du 

do. 



e i& (yn+Xyt>) 



(2.23a) 



(2.23b) 



Of course, different parametrization of the traveling-wave so- 
lution will also result in different combinations of modes. The 
shape of the modes in Eqs. (2.21) is shown in Fig. 1 for 5 = 
(NLS) and s = 2. 

In order to quantify the effect of perturbations, it is necessary 
to also employ the adjoint modes. Introducing the inner prod- 
uct as 

(y,w) =Refy*(x)w(x)dx= [(y Re w Re +y lm w lm )dx, (2.24) 
the adjoint of A[y, U] is found to be 



A^[y,U} 



2 d£ 2 



-A*y 



- 2i n u ^H") u kH'+^") y ^H') R ^'X") d %' d %" 

+ i JJ U {$+$>) U ($+Z")y($+Z>+Z") R (%'4") d % d % ■ (2-25) 

Note that A^ [y,U] = —iA[iy,U] . Using this property, one can 
immediately obtain the adjoint modes of Eqs. (2.21), i.e., the 
eigenmodes of A^: 



.du_ 



-2s- 



dU 



y x = \t>U, (2.26a) 



u, 



(2.26b) 



ds 
i dU 
A ~dt[ 

satisfying relations analogous to Eqs. (2.22). (Note, however, 
that ^ ±iyQ> an d the adjoint of a neutral eigenmode is a 
generalized eigenmode, and viceversa.) Due to the real and 
even nature of f(x) , the modes (2.21) and their adjoints (2.26) 
form a bi-orthogonal set which provides a basis of the gener- 
alized nullspace of A. Explicitly, 



where 



2s§-) E /A, 



(y A ,yA) = (y# 

(y x ,yx) = (y a ,yn) = \e/a, 



(2.27) 

(2.28a) 
(2.28b) 



and where E = \\u\\ 2 = J \ u(x,t)\ 2 dx is the pulse energy. Note 
that, unlike the case of the NLS equation, these adjoint modes 
are not normalized, and in general the norms ||ygl| 2 — (y^i^s) 
depend on both s and X, as shown in Fig 2. Of course the 
orthonormality could be achieved by properly redefining the 
adjoint modes, but the present choice is convenient for our 
purposes. For the NLS equation it is simply s = and E = 2A, 
so in that case the modes are indeed bi-orthonormal. As we 
show next, these linear modes and adjoint modes can be used 
to quantify perturbation-induced solution parameter changes. 
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FIG. 2: Top: the L2-norm of of the amplitude mode and the in- 
ner product (y^^A) as a function of the amplitude parameter A for 
different values of map strength. Bottom: the L2-norm of and the 
inner product (y^,y^>) as a function of the amplitude parameter A. 

C. Perturbation-induced parameter changes 

We now consider perturbations which manifest as an addi- 
tional term in the NLS equation (2.1): 



du 
dt 



dhi 
] d7 



'37 + M f A«)33 +g(t/t a )\u\ z u = ieh(x,t) , (2.29) 



with < £ <C 1 . Using the same multiple scale analysis as for 
Eq. (2.1), from Eq. (2.29) one obtains a perturbed DMNLS 
equation in which an inhomogeneous term is added to the 
right-hand side (RHS): 

du , -d 2 u 

+ II u (x+x') u {x+x") u *(x+x'+x") R {x',x") dxdx = ish(x,t) . 

(2.30) 

This formulation is of course very general, and it includes 
most physically interesting situations such as damping, am- 
plification, third-order dispersion, shock and Raman effects 
(e.g., see Refs. [6, 13]). 



Suppose u £ = u + ev solves Eq. (2.30), where u is a traveling- 
wave solution of the unperturbed DMNLS equation given by 
Eq. (2.12) and ev is the perturbation-induced solution change. 
Then the perturbation v(x,f) satisfies the perturbed linearized 
DMNLS equation, 



L[v] = h . 



(2.31) 



Of course, in order for the perturbation expansion to re- 
main well-ordered, the solutions of Eq. (2.31) must remain 
bounded. If, however, the right-hand-side of Eq. (2.31) has a 
non-zero component along the nullspace of L, secular terms 
will arise. As usual, such terms are removed by requiring that 
the parameters of the unperturbed solution become slowly- 
varying functions of time. Namely, introducing the fast and 
slow time scales t \ = t and ?2 = £t, one obtains 



du 
~dt 



du 
dti 



y-, du dQ 



(2.32) 



where now Q — A,£1,X,4>. The perturbed linearized DMNLS 
equation (2.31) therefore becomes 



j r l i V" du d Q i 



(2.33) 



where L\[v] is given by Eq. (2.14) with djdt replaced by 
d/dt\. Recalling Eqs. (2.19) and (2.23), one can rewrite 
Eq. (2.33) as 



dw . . dQ. 

drr Al[w]+Xy ^ 



Q " l 2 



(2.34a) 



where w(x,t) = e~' & ( x,t K(x,t) and where Ai[w] is given by 
Eq. (2.20) with t 1 = t\ and ^ = x — X (t \ ) . [The extra term pro- 
portional to y$ comes from Eq. (2.23b).] We now decompose 
the right-hand side of Eq. (2.34a) as 



Q 



(e m y 



Owe) 



so that 



y<2 



s ) = o 



-hn 



(2.35a) 



(2.35b) 



for Q — A,Q.,X ,<P. The solvability condition (namely the 
requirement that the resulting PDE for w contain no secular 
terms) then provides evolution equations for the solution pa- 
rameters: 



dA 
dt 

dQ. 

dt 



= e 



(y„,yn) 



dX 
~dt 



h) 



(2.36a) 



(2.36b) 



(2.36c) 



dt 



(A 2 -H 2 ) + e 



eX 



(e i& y_ a ,h) 
(la^yn) 



(2.36d) 
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In the special case h(x,t) = Au(x)8(t — t a ), Eq. (2.30) de- 
scribes changes in the initial condition. [Integrating from 
t = t — At to t = t + At and letting At — > one has u{x,t„) = 
u(x,tg)+eAu(x).] In this case, Eqs. (2.36) yield the parame- 
ter changes as Q(t+) = Q{t~) + £ AQ, where 



Owe) 



fovQ=A,£l,X, while 
A* = 



-x 



(e i@ y_ a ,Au) 

(ya,ya) 



(2.37a) 



(2.37b) 



All of these results reduce to the standard perturbation theory 
for the NLS equation (e.g., see Refs. [14, 24]) when s = 0. The 
connection between invariances, conservation laws, and linear 
modes of the DMNLS equation is further explored in Ref. [9]. 
In section III we will use the linear modes and Eq. (2.37) to 
guide the simulations of noise-perturbed lightwave systems. 



III. RARE EVENTS IN DISPERSION-MANAGED 
LIGHTWAVE SYSTEMS 



As mentioned in the introduction, the performance of both op- 
tical communication systems and certain femtosecond lasers 
is limited by noise. As a special but physically important 
example, we therefore now use the results of section II to 
quantify the effects of noise on dispersion-managed solitons. 
Specifically, we take h(x,t) in Eqs. (2.29) and (2.30) to be 



h{x,t) = £ v„(x)8(t-nt a ), 



(3.1) 



n=\ 



where t a is the dispersion map period, as before, and 8(t) is 
the Dirac delta distribution. In other words, u(x,t) evolves 
according to the unperturbed DMNLS equation (2.5a) at all 
times except at t = nt a , for n = l,...,N a , when 



u(x, rit^) = M ( x > nt a ) + e v « ( x ) > 



(3.2) 



where v„ (x) is taken to be normalized white Gaussian noise, 
satisfying 



E[v„(jc)] - 0, E[v„(x)v„*,(x')] = 8{x-x') 8 m 



(3.3) 



where E[-] denotes ensemble average, and where in this case 
the small parameter e 2 is the dimensionless noise variance. 
Starting from a soliton input pulse, namely u(x,0) given by 
Eq. (2.12), we are then interested in computing the probability 
density function (pdf) of the soliton parameters at the output 
time N a t a - 

Even if the statistics of the noise sources are Gaussian, the 
resulting statistics at the output is not, in general, because 
propagation is nonlinear. Indeed, as mentioned earlier, the 
combination of noise and nonlinearity presents a formidable 
challenge if one is interested in calculating the probabilities 
of errors when performance standards dictate that errors be 



rare events. It has recently been shown that variance reduc- 
tion methods such as importance sampling can be used to cal- 
culate pdfs in such systems accurate to very small probabili- 
ties [10, 23, 25, 35]. Here we use the results of section II to 
implement IS methods for the DMNLS equation in the pres- 
ence of noise. 



A. Importance sampling for the DMNLS equation 

The idea behind importance sampling is straightforward: in 
order to calculate the probability of a desired rare event, sam- 
ple the noise from a biased distribution that makes the rare 
events occur more frequently than would naturally be the case, 
while simultaneously correcting for the biasing. 

Consider a set of random variables X = (jq, . . . ,xn) dis- 
tributed according to a joint probability distribution p(X). The 
probability P that a function y(X) falls into a desired range 
can be expressed via the multi-dimensional integral 

P = P[>' G Yd] = E[/(y(X))] = fl(y(x))p(x)dx, (3.4a) 

where the indicator function I(y) equals 1 when y e Yj and 
otherwise. An unbiased estimator for P can be constructed 
via Monte-Carlo (MC) quadrature as 



M 



(3.4b) 



m=l 



where the M samples X m are drawn from the distribution 
p(X). If P is very small, however, an unreasonable number of 
samples are necessary to produce events for which y is in Lj, 
let alone enough to accurately approximate the integral. How- 
ever, one can rewrite Eqs. (3.4) as follows: 



F\yeY d ] = fI(y(x))L(x)p*(x)dx, 

j M 

m= 1 



(3.5a) 
(3.5b) 



where the samples X* i/n are now drawn from the biasing dis- 
tribution p*(X), and where the quantity L(X) = p(X) / p*(X) 
is the likelihood ratio. When an appropriate biasing distri- 
bution is selected, importance-sampled Monte-Carlo (ISMC) 
simulations can accurately estimate the probability of the 
sought-after rare events with a small fraction of the number 
of samples that would be necessary with straightforward MC 
methods. The challenge is, of course, to properly choose the 
biasing distribution. Indeed, in order for importance sampling 
to work, p* (X) should preferentially concentrate the MC sam- 
ples around the most likely system realizations that lead to the 
rare events of interest. In our case the random variables X are 
the noise components added at the end of each dispersion map 
period. Thus, in order to successfully apply IS we must find 
the most likely noise realizations that lead to a desired value 
of the soliton parameters at the output. 

By substituting V„(x) into Eq. (2.37) of section II C, we im- 
mediately obtain the noise-induced parameter change at the 
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n th map period as: 



Ref e - i& Wf Q {x)v n {x)dx 



(3.6) 



where Q — A,Q.,X as before, with a slightly more complicated 
expression for A4>„. Moreover, for white Gaussian noise, 
maximizing the probability of a specific noise realization is 
equivalent to minimizing the negative of the argument of the 
exponential in the pdf; that is, minimizing the quantity 



J\Vn( X )\ 2 d X . 



(3.7) 



Hence, in our case, the problem of determining the optimal bi- 
asing amounts to finding the noise realization that minimizes 
the integral in Eq. (3.7) subject to the constraint of achieving 
a desired parameter change, that is, subject to the constraint 
AQ„ = Agtarget, with AQ n given by Eq. (3.6). This optimiza- 
tion problem can be solved by formulating Eqs. (3.7) and (3.6) 
as a Lagrange multiplier problem, as in Refs. [23, 29]. Solv- 
ing the resulting problem then yields V„(x) — V„ i0 pt(*)> where 



ReJy* Q (x')y Q (x')dx' ^ 

v "' optW - J\y_ Q (x>)\idx> Aetargete 



y_ Q (x). (3.8) 



To induce a larger than normal parameter change, we then 
concentrate the MC samples around this optimal path. We 
do so by biasing the noise adding V„ j0p t(x) as a deterministic 
component; that is, we take 



V*,n( x ) = v «,optW + V n (x) . 



(3.9) 



where V„ ]0p t(^) is given by Eq. (3.8), and where V„(x) is also 
a white noise process satisfying Eqs. (3.3). 

Note that the optimal path V„ j0p t(*) depends on both the eigen- 
modes and the adjoint eigenmodes of the linearized DMNLS 
equation found in section II. In particular, Eq. (3.8) implies 
that the optimal biasing is proportional to the adjoint eigen- 
mode of the quantity that one desires to change, a result that 
might not be obvious apriori. 

Note also that, once the most likely noise realization that 
produces a given parameter change AQ n at each map period 
are known, one should also find the most likely way to dis- 
tribute a total parameter change AQ tot at the output among 
all map periods (cf. Refs. [9, 24, 29]). This further op- 
timization problem can also be formulated and solved as a 
variational problem [9]. When targeting large amplitude or 
frequency changes, however, it suffices to simply distribute 
equally this total change among all map periods. Thus, in 
the numerical simulations described in section IIIB we set 
Agtarget = AQ tol /N a for n = 1 , . . . ,N a . 



B. Importance-sampled Monte-Carlo simulations 

We now apply the ideas presented above to concrete numeri- 
cal experiments of dispersion-managed systems under the ef- 
fect of noise. We perform importance-sampled Monte-Carlo 



(ISMC) simulations of the DMNLS equation (2.5) perturbed 
by noise, and we compare the results to standard Monte-Carlo 
simulations of the original NLS equation with dispersion man- 
agement (2.1), also subject to noise. 

Let us discuss the approach we used for the numerical simu- 
lations of the noise-perturbed DMNLS equation. We simulate 
the evolution of a dispersion-managed optical signal by solv- 
ing Eq. (2.5b) numerically and adding noise to the signal at pe- 
riodically spaced time intervals. The initial condition, i.e., the 
input DMS shape at t = 0, is generated by solving the nonlin- 
ear integral equation Eq. (2. 10) as explained in Appendix A 2, 
and time evolution is performed with a fast numerical method, 
as described in Appendix A 1 . White noise, which is added 
to the signal at t = nt a for n = l,...,N a , is numerically dis- 
cretized as a collection of independent, identically distributed 
zero-mean normal random variables, one each for the real and 
imaginary parts of the signal at each spatial grid point. Propa- 
gation and the addition of noise continue in this way until the 
signal reaches the output at f out = N a t a . 

In standard Monte-Carlo simulations, one repeats the above 
process for several different noise realizations while monitor- 
ing the output value of the quantities of interest (e.g., energy 
and/or frequency), and then computes their statistics. 

For importance-sampled Monte-Carlo simulations, one also 
uses the basic framework described above. If one wants to 
obtain larger-than-normal deviations of a quantity Q, however, 
one also performs the following steps at each map period be- 
fore adding the noise: 

1 . Determine the underlying DMS from the noisy signal. 
We do this by using the noisy pulse as the initial con- 
dition in the iteration scheme for the nonlinear integral 
equation (A. 5), where the denominator in the RHS is 
replaced by X +d(co — Q.) 2 to automatically obtain a 
DMS shape with the correct carrier frequency. The val- 
ues of Q. and X are obtained from the (low-pass fil- 
tered) noisy signal, respectively by using moments and 
by matching the total energy in a pre-computed table of 
DMS energy as a function of X. 

2. Obtain the linear modes and adjoint modes of the lin- 
earized DMNLS equation around the given DMS. We 
do this by numerically calculating the x and s deriva- 
tives of the underlying DMS, and then using Eqs. (2.21) 
and (2.26). The x-derivative is calculated using pseudo- 
spectral methods, while the derivative with respect to s 
is calculated by performing step 1 twice, once at the 
given value of s and once at s + ds. 

3. Generate an unbiased noise realization, shift its mean 
with the appropriately scaled adjoint mode associated 
with Q according to Eqs. (3.9) and (3.8), and update the 
likelihood ratio. 

One then adds the noise to the pulse, propagates the noisy 
signal to the next map period, and repeats this process until 
the signal reaches the output. For a given simulation, sev- 
eral thousand ISMC samples, generated with a few different 
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FIG. 3: Samples from ISMC simulations of the DMNLS equation. 
Here, the pulse energy (normalized to input energy) is plotted as a 
function of time (i.e., distance in physical units). The arrows rep- 
resent the different targeted output energies: a larger than normal 
output energy (blue), a smaller than normal output energy (red), and 
unbiased energy (black). Also plotted are deterministic paths (thick, 
smooth curves, with color corresponding to the target) predicted by 
our perturbation theory. These are the preferential paths around 
which we attempt to sample by biasing the noise with the adjoint 
linear modes. For each of three different targeted output energies, a 
few dozen ISMC samples are also shown (also colored correspond- 
ingly), demonstrating that the actual trajectories indeed follow the 
predictions of the theory. See text for a detailed discussion of system 
parameters. 

biasing targets, are collected, and their contributions are com- 
bined using multiple importance sampling [10, 36] in order to 
numerically generate the pdf of the quantity of interest. 

Even though at each map period the noise induces only a small 
change in the solution parameters, these small changes can 
accumulate into large parameter changes at the output, result- 
ing in a significantly distorted signal. This gradual build up 
of noise-induced distortions is evident in Figure 3, where we 
plot the energy as a function of time for several different noise 
realizations biased around the optimal energy path predicted 
by the theory. For each sample path in this figure, the noise 
added at each map period is biased by adding a proper multi- 
ple of the adjoint amplitude mode in order to change the pulse 
amplitude and hence its energy. Note how the random samples 
are concentrated near the trajectories predicted by the theory. 
One can think of these trajectories as a low-dimensional pro- 
jection of a near-optimal path through state space to reach a 
targeted rare event. 

Note also that the linearized DMNLS equation is only used 
to guide the biasing, while, for each individual sample in 
the ISMC simulations, the full DMNLS equation is solved to 
propagate the signal. That is, no approximations are used for 
propagating the signal, and no assumptions are made about its 
shape to predict or calculate the pulse parameters at output. In 
other words, the only approximation in the simulations (be- 
yond roundoff and truncation due to discretization) lies in us- 
ing the information based on the linearized DMNLS equation 
in order to bias the noise. Thus, use of importance sampling 
enables full nonlinear simulation of large, noise-induced pa- 
rameter changes. 



C. Results and discussion 

As a test of the method, we performed numerical simulations 
of two noise-perturbed dispersion-managed systems with dif- 
ferent sets of physical parameters. For both systems we com- 
pared ISMC simulations of the DMNLS equation (2.5) to 
standard MC simulations of the original NLS equation with 
dispersion management (2.1), numerically integrated with a 
standard second-order split-step Fourier method. For both 
systems we take a piecewise constant dispersion map. That 
is, we consider the transmission link to be comprised of alter- 
nating sections of fiber with opposite signs of dispersion. 

In the first system, which we will refer to as system (a), we 
seek large changes in frequency at the output. In the sec- 
ond system, referred to hereafter as system (b), we seek large 
changes in amplitude and hence in output energy. In both 
cases we choose the system parameters based on realistic val- 
ues for optical fiber communication systems. Typical values 
of system parameters for the DMNLS as a model of fem- 
tosecond lasers can be obtained from Ref. [28]. We consider 
an average dispersion of 0.15ps 2 /km a nonlinear coefficient 
of 1.7 (W-km)- 1 , and a fiber loss of 0.21 dB/km. We set a 
unit time of 17 ps, normalizing x in Eq. (2. 1) correspondingly, 
and we use the resulting dispersion length of 1,923 km to nor- 
malize t, resulting in d = 1 in Eq. (2.1). We further space 
amplifiers every 100 km (resulting in t a = 0.052), taking the 
dispersion map period to be aligned with the amplification 
period. The corresponding power needed to have g = 1 in 
Eq. (2. 1) is 2.96 mW, and we use this value as a unit to normal- 
ize pulse powers. For system (a) we consider a map strength 
of s = 2 and a propagation distance of 2,000 km (resulting 
in N a = 20). For system (b) we consider a map strength of 
5 = 4 and a propagation distance of 4,000 km (resulting in 
N a = 40). Thus, in both systems the average dispersion is 
small, while the local dispersion is large in magnitude. [Re- 
call that the map strength parameter s quantifies the difference 
in magnitude between local and average values of dispersion, 
cf. Eq. (2.7).] For system (a) we take X — 1.5 (correspond- 
ing to an initial pulse with 6.66 mW peak power) and for sys- 
tem (b) we take X = 2 (corresponding to 11.8mW). Finally, 
assuming a spontaneous emission factor of 1.5, system (a) has 
an optical signal-to-noise ratio of 16.7 (resulting in a dimen- 
sionless noise variance £ 2 = 2.372 • 10~ 4 ), system (b) of 13.8 
(resulting in e 2 = 9.486 • 10~ 4 ). 

The output distributions of both frequency and energy are of 
course of practical interest in communications, since large de- 
viations of each quantity will result in transmission errors. 
Frequency changes translate in group velocity changes, and 
hence in the pulse walking off its assigned bit slot. Similarly, 
a pulse that loses a significant fraction of its energy will be 
incorrectly detected in an amplitude-shift-keyed system. 

In Figure 4 we plot the pdf of the output frequency of pulses 
from the dispersion-managed system (a) described above. We 
perform standard MC simulations of the noise-perturbed NLS 
equation with DM (2.29), with 100,000 samples to estimate 
the pdf of output frequency (red dots). We also plot a Gaus- 
sian pdf (black dashed line) whose variance is consistent with 
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FIG. 4: Probability density function of output frequency in a 
dispersion-managed system. The solid cyan curve shows results from 
ISMC simulations of the DMNLS Eq. 2.5 with 75,000 samples. The 
red dots result from standard MC simulations of the NLS equation 
with DM, Eq. (2.1) with 100,000 samples. The black curve is a 
Gaussian pdf obtained from Gordon-Haus theory of the NLS equa- 
tion with DM. Note how unbiased MC simulations of the NLS equa- 
tion with DM agree exactly with the ISMC simulation of the DMNLS 
and the Gaussian fit to that simulation as far down in probability as 
the unbiased simulations can reach. 



solid line). The results of ISMC simulations of the DMNLS 
equation, of standard MC simulations of the NLS equation 
with DM, and the Gaussian fit to that simulation all match 
exactly to very small probabilities. This comparison demon- 
strates the effectiveness of using the modes of the linearized 
DMNLS equation to find rare events in dispersion-managed 
optical systems. 

In Figure 5 we plot the pdf of the output energy of pulses in 
dispersion-managed system (b) described above. Here, as in 
Fig. 3, energy is normalized by the the energy of the input sig- 
nal, i.e., by the "back-to-back" signal energy. Again, the red 
dots show results of standard MC numerical simulations of 
the noise-perturbed NLS equation with DM, with 1,000,000 
samples; the black dashed line shows a Gaussian fit to the re- 
sults of these simulations, and the cyan solid line shows the 
results of ISMC simulations of the noise-perturbed DMNLS 
equation, with 42,000 samples. It is worth noting that the pdf 
generated from standard MC simulations of the NLS equation 
with DM clearly deviates from the Gaussian fit, but it agrees 
very well with the ISMC simulations of the DMNLS equation, 
as far down in probability as the unbiased MC simulations 
can reach. These comparison provides a strong validation of 
the DMNLS equation as a model of dispersion-managed light- 
wave systems. 
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FIG. 5: Probability density function of normalized output energy. 
The solid cyan curve shows results from ISMC simulations of the 
DMNLS with 42,000 samples. The red dots result from standard 
MC simulations of the NLS equation with DM with 1,000,000 sam- 
ples. The black curve is a Gaussian fit to that simulation. Note how 
unbiased MC simulations of the NLS equation with DM clearly de- 
viate from Gaussian, but agree well with ISMC simulations of the 
DMNLS as far as down in probability as the unbiased simulations 
can accurately reach. 

a theoretical model of Gordon-Haus effect for a DM system 
governed by the noise-perturbed NLS equation in the pres- 
ence of dispersion management (2.29) [16, 21, 32], assuming 
a Gaussian ansatz for the pulse shape at the chirp-free point. 
Finally, we perform ISMC simulations of the noise-perturbed 
DMNLS equation (2.30) with 75,000 samples, and we plot the 
the corresponding results for the pdf of output frequency (cyan 



IV. CONCLUSIONS 

We described a perturbation theory for soliton-based dis- 
persion-managed lightwave systems, whose dynamics is gov- 
erned by the dispersion-managed NLS (DMNLS) equation, 
and we used the results of the theory to guide importance- 
sampled Monte-Carlo simulations to quantify the effects of 
noise in these systems. The present theory differs from the 
soliton perturbation theory which applies to the NLS equa- 
tion in several important respects. First, due to the loss of 
integrability, the eigenmodes of the linearized DMNLS equa- 
tion are derived from the invariances of the equation rather 
than from the inverse scattering method. Second, the DMNLS 
equation is not scale-invariant, but is invariant under a gen- 
eralized scaling transformation; as a consequence the ampli- 
tude mode depends explicitly on the map strength. Third, un- 
like for the NLS equation, the linear modes of the DMNLS 
equation are not automatically bi-orthonormal, and in partic- 
ular their norms and inner products depend on both the map 
strength and the pulse energy. 

The results of importance-sampled numerical simulations of 
the noise-perturbed DMNLS equations agree very well with 
the results of Gordon-Haus theory for dispersion-managed 
systems (which are based on the original NLS equation) as 
well with the results of standard Monte-Carlo simulations of 
the original NLS equation with dispersion management as far 
down as those can go in probability. This is true even when 
those results deviate from Gaussian distributions. Both of 
these results provide a further important test of the validity 
of the DMNLS equation as a model of dispersion-managed 
lightwave systems. 
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It should be noted that, in some parameter regimes, the 
DMNLS equation also admits some internal modes [11, 26]. 
The accumulation of noise components in such internal modes 
could also contribute to system failures. If no generalized 
modes are associated with these modes, however, the vari- 
ance of the resulting noise-induced pulse fluctuations would 
grow linearly in time, as opposed to cubically (as is the case 
for phase and timing fluctuations). Therefore, one would ex- 
pect these modes not to be the dominant source of errors. Of 
course, this argument should be validated by more precise cal- 
culations and/or careful numerical experiments. We will ad- 
dress this issue in the future. 
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APPENDIX: FAST NUMERICAL METHODS FOR THE 
DMNLS EQUATION 

Here we discuss fast numerical methods for the DMNLS 
equation (2.5) which were used in our numerical simulations. 
We address two issues: (i) numerical methods for time in- 
tegration, and (ii) numerical methods to find traveling-wave 
solutions. 



1. Time integration of the DMNLS equation 

We first discuss time integration methods for the DMNLS 
equation. Equation (2.5a) differs from a standard PDE be- 
cause of the double integral, and it should be obvious that the 
most computationally expensive task when trying to integrate 
it numerically is the evaluation of the double integral in either 
of Eqs. (2.5). If N x is the number of grid points in x (or a>), 
a straightforward quadrature scheme requires 0(N X ) opera- 
tions (since one needs to evaluate a different double integral 
for each value of x or a>). As we now show, however, it is 
possible to evaluate the integral with only 0(JN x \ogN x ) op- 
erations, where J is in principle an arbitrary constant whose 
meaning will become clear shortly. 

Let us denote the double integral in Eq. (2.5b) as 

K(co,t) = X^M( ffl+ra0 M( ffl+(a » ) M(; +(i) / +(i) « ) '-( (a / £l) « ) da)'da)' , 

(A.l) 

(where we reinstated the primes to distinguish the solution of 
the DMNLS equation from that of the original NLS problem). 
Recall from Eq. (2.3) that u((0,t,Q = u'(co 1 t)e- iC( ^°> 2 / 2 
is the leading-order solution of the original NLS equation 
with dispersion management, namely Eq. (2.1). Also recall 



from Eq. (2.6a) that the kernel r(x) is an average over a dis- 
persion map period, and R(t,t') = ^""/^((Ofi)')]. Indeed, re- 
tracing backwards the steps of the multiple scale expansion 
used to obtained the DMNLS equation [1], one can express 
the whole quantity K(co,z) as an average: 

K(co,t) = )e c ^ a2 l 2 g{Q ri0 [(\u\ 2 u)(x,t, Q] dC , (A.2) 
o 

where u(t 7 z 7 Q = 7t~ l [&((D,z, £)]. Note that t is a constant in 
the integral on the RHS of Eq. (A.2), and this is precisely the 
key for the fast numerical computation of K((o,t). 

Divide the interval [0, 1] into J equally spaced points 
£o, • • • , £/> w i tn Co = and Q = 1. For each fixed value of t 
we can calculate K(co,t) as follows: 

1. Fix £j and construct u(co,t, £ 7 ) = u' ((0,t)£~ lC ^ a ' 2 1 2 

2. Take the inverse Fourier transform to obtain u(x,t,£j) 
and construct the product \u\ 2 u. 

3. Take the direct Fourier transform to obtain iFa>[|w| 2 w]. 

4. Multiply the result by g(£ / )e' c( £' )£o2/2 to obtain the in- 
tegrand in Eq. (A.2). 

5. Repeat steps 1-4 for all grid points and evaluate the 
integral in Eq. (A.2) to obtain K((0,t). 

The DMNLS equation can now be integrated in time using any 
desired numerical scheme. For example, one can use an exact 
integrating factor on the linear part of Eq. (2.5) and an explicit 
fourth-order RK4 for the nonlinear part. A few remarks are 
now in order: 

o The above procedure, which can be carried out for any 
choice of dispersion map and nonlinear coefficient, is a 
generalization of an algorithm originally introduced in 
Ref. [20]. 

o Steps 2 and 3 each cost 0(N x logN x ) operations. Since they 
must be repeated for each of the grid points £o, • • • , £/, the 
overall complexity of the algorithm is 0(JN x logN x ). 

o In essence, the algorithm reduces the calculation of K(co,z) 
to that of the effective nonlinearity in the original NLS 
equation, averaged over one period of the dispersion map. 

o In practice, the value of J is dictated by the value of 
the map strength and the need to adequately resolve the 
changes in the integrated dispersion function C(Q. The 
same requirements however also dictates the integration 
step size in the original NLS problem. The computational 
complexity of the DMNLS equation (2.5) is thus the same 
as that of the original NLS problem (2.1). 

o It appears more advantageous to integrate Eq. (2.5b) rather 
than Eq. (2.5a), since this allows one to treat the linear 
(stiff) portion of the PDE exactly. Also, the use of a split- 
step method does not seem as desirable here (unlike the 
case of the NLS equation), since it is not possible to inte- 
grate the nonlinear portion exactly. 
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2. Traveling-wave solutions of the DMNLS equation 

We now discuss a numerical method to find stationary solu- 
tions of the DMNLS equation, i.e., dispersion-managed soli- 
tons (DMS). That is, we look for solutions of Eq. (2.5a) in the 
form 

t t A {x,t)=f{x)e &h l 2 , (A3a) 
or, equivalently, solutions of Eq. (2.5b) in the form 

*4(©,0=/(o))e a2 ' /2 . (A.3b) 



Once such a stationary solution is available, a three-parameter 
family of travelling waves can be generated using the transla- 
tion, phase and Galilean invariance of the DMNLS equation. 
That is, if u' st (x,t ) is any solution of Eq. (2.5a), so is 

^ ^t) = ^- da2t l 2+ ^u' sl {x-x -d£lt,t), (A.4) 

where £2, x and <p are arbitrary real parameters. 

Inserting Eq. (A. 3b) into the DMNLS equation yields the non- 
linear integral equation (2.10), which we rewrite here for con- 
venience: 

/(«) = 

2 - 

X 2 + dco 2 ■ttf( m + m 'rf( m + m '')f*®+<»'+(>>'') r ( co ' a '') d(0 ' d(0 " ■ 

(A.5) 

For each fixed value of X, Eq. (A.5) yields the shape of the 
corresponding dispersion-managed soliton. (Or, rather, its 
Fourier transform.) Thus, for each value of s, X plays the role 
of a nonlinear eigenvalue, which is in one-to-one correspon- 
dence with the dispersion-managed soliton's energy. For the 
NLS equation, one simply has u(x,t) — Asech[Ax]e iA 'I 2 ; thus 
X = A is exactly half of the pulse energy: / | u s (x, t ) 2 dx = 2A . 
This is related to the existence of a simple scaling invariance: 
if u a (x,t) is any solution of NLS, so is u(x,t) = X u (Xx,X 2 t). 
When the map strength is non-zero, however, this invariance 
is lost, and the scaling invariance of the DMNLS equation 
is more complicated than that of the NLS equation, as dis- 
cussed in section II. As a consequence, a different integral 
equation (A.5) must be solved to obtain the soliton shape for 
each given value of X, unlike the NLS equation. 

Let R[f(a>)] denote the right-hand side of Eq. (A.5). A naive 
approach to solving Eq. (A.5) is to simply apply a Neumann 



iteration scheme: — ^[/S]- This iteration scheme is 

divergent, however. The key is to apply a modified iteration 
scheme: 



(A.6) 



where the convergence factor C is used: 

S L [/(CO)] 



C[f {m) ] 



«[/(0 



(A.7a) 



with 



s l\f{a>)\ = / l/(co)l dC0, SR\f(a>)\ = //(*»)#[/(<»)] dtO 



(A.7b) 



Again, a few remarks are in order: 



o The method discussed in section A 1 for the fast numerical 
computation of the double integral in the DMNLS equation 
also applies, of course, for the calculation of the double 
integral in Eq. (A.5). 

o Since C[/( ffl )] = 1 for all solutions of Eq. (A.5), any 
solution of Eq. (A.5) is also a solution of the modified inte- 
gral equation f {m) = \C[f {m) \ \ a R[f( m) \, for which Eq. (A.6) 
is a standard Neumann iteration scheme. 

o Since C[f^\ — 1 when solves Eq. (A.5), the value 
of C[/( ffl )] can be used as a monitor of convergence, for 
example by requiring that 1 1 — C[/( ffl )] | drop below a pre- 



defined threshold (e.g., 10 or 10 ) as a termination 
condition. 

o This iteration scheme, which was first used in Ref. [1] to 
find the shape of dispersion-managed solitons, is based on 
a method introduced in Ref. [27]. A proof of the conver- 
gence of this method for a general class of nonlinear evo- 
lution equations was recently given in Ref. [31]. 

o As might be expected, the choice of a is crucial. Indeed, it 
can be shown that the method converges for 1 < a < 2, and 
optimal convergence is obtained for a = \ . Note, however, 
that d > is required for convergence. If d < instead, the 
denominator of the RHS of Eq. (A.5) has two simple poles 
at ft) = ±X I \[\&\, and even the modified iteration scheme 
diverges. 
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